/* +++Date last modified: 05-Jul-1997 */

/*
**  CUBIC.C - Solve a cubic polynomial
**  public domain by Ross Cottrell
*/

#include <math.h>
#include <stdlib.h>
#include "snipmath.h"

void SolveCubic(double  a,
                double  b,
                double  c,
                double  d,
                int    *solutions,
                double *x)
{

      double a1 = d1_div(b, a);
      double a2 = d1_div(c, a);
      double a3 = d1_div(d, a);
  //long double    a1 = b/a, a2 = c/a, a3 = d/a;
      double l1a = d1_mul(a1, a1);
      double l1b = d1_mul(-3.0, a2);
      double l1c = d1_add(l1a, l1b);
      double Q = d1_div(l1c, 9.0);
 //long double    Q = (a1*a1 - 3.0*a2)/9.0;
      double l2a = d1_mul(2.0, a1);
      double l2b = d1_mul(l2a, a1);
      double l2h = d1_mul(l2b, a1);
      double l2c = d1_mul(-9.0, a1);
      double l2d = d1_mul(l2c, a2);
      double l2e = d1_mul(27.0, a3);
      double l2f = d1_add(l2h, l2d);
      double l2g = d1_add(l2f, l2e);
      double R = d1_div(l2g, 54.0);
 //long double R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0;
      double R2 = d1_mul(R, R);
      double Q3a = d1_mul(-Q, Q);
      double Q3 = d1_mul(Q3a, Q);
      double R2_Q3 = d1_add(R2, Q3);
//double    R2_Q3 = R*R - Q*Q*Q;

      double    theta;

      if (R2_Q3 <= 0)
      {
            *solutions = 3;
            theta = acos(R/d_sqrt(Q*Q*Q));
            x[0] = -2.0*d_sqrt(Q)*cos(theta/3.0) - a1/3.0;
            x[1] = -2.0*d_sqrt(Q)*cos((theta+2.0*PI)/3.0) - a1/3.0;
            x[2] = -2.0*d_sqrt(Q)*cos((theta+4.0*PI)/3.0) - a1/3.0;
      }
      else
      {
            *solutions = 1;
            x[0] = d_pow(d1_add(d_sqrt(R2_Q3),fabs(R)), d1_div(1,3.0));
            x[0] += Q/x[0];
            x[0] *= (R < 0.0) ? 1 : -1;
            x[0] -= a1/3.0;
      }
}

#ifdef TEST

int main(void)
{
      double  a1 = 1.0, b1 = -10.5, c1 = 32.0, d1 = -30.0;
      double  a2 = 1.0, b2 = -4.5, c2 = 17.0, d2 = -30.0;
      double  x[3];
      int     solutions;

      SolveCubic(a1, b1, c1, d1, &solutions, x);

      /* should get 3 solutions: 2, 6 & 2.5   */

      SolveCubic(a2, b2, c2, d2, &solutions, x);

      /* should get 1 solution: 2.5           */

      return 0;
}

#endif /* TEST */
